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Abstract. The problem of determining a metal’s Fermi surface from measurements of 
projections of the electron or electron/positron momentum densities, such as obtained 
by Compton Scattering or Angular Correlation of Positron Annihilation Radiation, 
respectively, is considered in a Bayesian formulation. A consistent approach is presented 
and its advantages compared to previous practice is discussed. A validation of the 
proposed method on simulated data shows its systematic accuracy to be very satisfactory 
and its statistical precision on modest experimental data to be surprisingly good. 
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1. Introduction 

The classical approach for experimentally determining the Fermi surface of metallic 
systems is to exploit its effect on quantum oscillations, such as in the de Haas-van Alphen- 
or the Shubnikov-de Haas-effects. These techniques, while proven to be very powerful 
for obtaining quantitative information on the dimensions of the Fermi surface with high 
precision, vitally depend on long scattering lengths of the electrons, and therefore are 
applicable only at cryogenic temperatures and at vanishing occupational disorder. Apart 
from that, the task of assigning the measured extremal orbits to actual features of the 
Fermi surface becomes challenging for systems with more complicated multi-sheet Fermi 
surfaces (e.g., Brasse et al. 2013). 

Experimental methods that measure electron momentum densities work equally well 
for disordered states (both occupational disorder and temperature) and therefore hold 
the promise of determining the Fermi surface from the discontinuities in the occupied 
momentum densities as long as the concept of a well-defined Fermi surface is meaningful 
at all (Dugdale 2014). These methods comprise Compton scattering (Cooper 1985) and 
Angular Correlation of Positron Annihilation Radiation (ACPAR, also ACAR, Bisson 
et al. 1982). However, in these techniques the primary experimental data are plane 
projections (Compton scattering, early positron annihilation) or line projections (recent 
positron annihilation setups), from which the three-dimensional momentum density has to 
be computationally determined. Also, in positron annihilation experiments the sampled 
two-photon momentum density differs from the underlying ideal electron momentum 
density due to positron wave-function effects and electron-positron correlations, although 
the position of the discontinuity due to the Fermi surface will remain unaffected. 
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The problem of reconstructing densities from projections has seen much attention 
due to its relevance for medical and technological imaging. For used approaches in the 
specific case of momentum densities in solids see, e.g., the recent review by Kontrym- 
Sznajd (2009). In short, densities were reconstructed in the majority of previous works 
either by discretizing an analytical inverse of the Radon transform (Radon 1917), or 
by expanding the measured projections and sought densities into basis functions with 
convenient transformation behaviour. Both of these approaches are direct methods 
in the sense that the results are derived by applying a sequence of explicitly defined 
transformations to the data. As a consequence, their computational complexity is modest, 
which historically was a reason for their adoption, and for diagnostical applications 
still is. On the other hand, experimental methods for determining electron momentum 
densities are countrate-limited, therefore it should be the power of an analysis method 
rather than the runtime which dictates the method to be preferred. 

In this paper, we will give a Bayesian formulation of the data analysis problem. For 
illustration, we will concentrate on the case of two-dimensional ACPAR, although our 
approach is equally applicable to Compton scattering. We will show that the formulation 
corresponds to a regularized inverse problem, and we will illustrate how its solution can be 
practicably obtained. The main features of our method lie in the avoidance of systematic 
errors thanks to a consistent description of the whole problem, and quantitative results 
due to an explicit parametrization of the Fermi surface. We will demonstrate the power 
of our method by applying it to simulated data, which will allow us to conclude that 
an ACPAR experiment with moderate statistics and resolution is able to determine 
the shape of the Fermi surface quantitatively with an accuracy that is comparable to 
quantum-oscillatory methods, at elevated temperatures and in the presence of disorder. 

2. Definition of the problem 

2.1. Bayesian formulation 

In the problem at hand, the basic unknown quantity is the three-dimensional momentum 
density p(p). For Compton scattering, this concerns the actual electron momentum 
density (Fourier components of the electron states), while in ACPAR the positron wave- 
function and electron-positron correlation effects modify the probed density, which is then 
termed two-photon momentum density. In either case, the density is a smooth function, 
apart from steps when crossing sheets of the Fermi surface. In the extended-zone scheme, 
the density will decay towards high p and display the point group symmetry of the 
crystal, while the Fermi surface features have the crystal’s space group periodicity. 

The experimental data y are given by a set of one- or two-dimensional spectra, 
corresponding to plane or line projections of the underlying three-dimensional densities. 
To a very good approximation, the data follow Poissonian statistics, in particular the 
noise for distinct data points is independent, and the noise probability distribution is 
uniquely defined by the expectation (and can be estimated from the measured signal). 
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The relation between these two quantities is given by a linear operator P, which 
specifies how the expected value of the measured signal (i.e., before Poisson quantization) 
results from a given, but experimentally a priori unknown, momentum density. We will 
call it the projection operator, as its function is to essentially integrate over the transversal 
or longitudinal momentum components for the distinct experimental orientations of 
the crystal. Additionally, it also takes the detector efficiency and momentum sampling 
functions and the smearing of the spectra due to finite resolution into account, which 
will be discussed in more detail in Sect. 13.21 

Thus, our initial formulation for the problem to be solved is the following: What is 
the posterior probability distribution p post for p, given an outcome of the experiment 
y and considering prior knowledge or assumptions, quantified in the prior distribution 
Pprior(p)? By Bayes’ formula, the answer is formally given by 

Ppost(ply) = Plit(gl ^p ( A (i) 

where puk(y | p) is the so-called likelihood function, the modelled probability for observing 
the actual experimental outcome y in a repeated experiment under the assumption of 
the momentum density p. In the present case it is just 

Phk(y\p) = ]^[ PPoisson \yii (P(p))jJ =: -^Poisson (lb P (p) ) j (2) 

i 

where the measured spectra are treated as a vector (y)* and 

PPoisson(^; A) e , (3) 

is the familiar expression for the Poisson distribution, with Pp 0 i SS on its version for vector 
arguments. 

The challenge to the physicist lies in formulating an expression for p pr io r (p) that 
takes into account the available understanding of the problem and therefore ensures a 
physically meaningful result. A first step towards this goal is the strict requirement for 
p to conform to the point symmetry group of the crystal. Instead of encoding this via a 
5-like part in p pr i or (p ) 5 if is more efficient to describe the densities only by their values in 
the irreducible wedge according to the point symmetry, with appropriate continuation 
over all of the momentum space. From here on, p is to be understood in this sense. 

The defining idea of our proposed approach actually follows just from letting 
modelling be guided by the physical picture: The expected behaviour of the density 
(smooth variations apart from jumps at the Fermi surface) is due to its being composed 
of contributions from the distinct conduction bands, each multiplied by a Fermi-Dirac 
occupation function with a practically discontinuous jump, sitting atop the contribution 
from the core states. The bare band densities (i.e., before taking occupation into account) 
will in fact be smooth, as can be derived from any simply tractable band structure 
model (be it nearly free electrons, tight binding, or the effective potentials in density 
functional theory). Instead of considering the total occupied density displaying jumps at 
the Fermi surface as the fundamental unknown, the most natural formulation therefore 
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is to have as free parameters both (1) a parametrization of the Fermi surface sheets in 
each conduction band (and spin channel in the case of magnetic ordering) and (2) the 
bare densities of the respective conduction bands plus the summed contribution from 
the core bands. The prior distribution p pr i 0 r(p) can then be used to favour smooth band 
densities with, e.g., additional positivity constraints or assumptions on the decay with p, 
or also information on the Fermi surface shapes, if available from prior experiments or 
calculations. Note that also for choosing the parametrization of the Fermi surfaces the 
physical picture indicates the natural way as the level set of a smooth function with the 
space group symmetry (in other words, the band dispersions). 

With these observations, we can rewrite Eq. (JT|) as 

Ppost(p,(r\y) oc P Poisson (y| Acr(p))p pr i 0 r(p, cr) with PX ff , (4) 

where p is now to be understood as the bare band densities and the action of X CT for given 
Fermi surface parameters a is to multiply the bare band densities by the occupations 
and sum over the band index, i.e., it is again a linear operator. 

If the primary quantity to be determined is the Fermi surface, Eq. Q can be 
marginalized over p to give the posterior distribution of the Fermi surface parameters 

PpostHy) = f dpp post (p, cr\y)- (5) 

Note that as we will elaborate below, the width of p pos t{p, cr| y) with respect to p does 
not vary much as a function of cr, so for practical purposes p pos t{v\y) is proportional to 
the maximum of p pos t(p, cr\y) for given a. 

2.2. Comparison to previous approaches 

The main points that distinguish our approach from the majority of those used previously 
are the following: first, it is formulated as a general problem of Bayesian inference instead 
of as an recipe of transforms obtained from analytical manipulations, and second, the 
Fermi surface is treated explicitly during the reconstruction instead of determined 
afterwards from the reconstructed densities. Here we will discuss the implications of 
these differences. 

The Radon transform is a bijection (in particular, it is invertible) between suitably 
regular n-dimensional function spaces (Radon 1917). Specifically, it corresponds to 
the relation between a pointwise defined function and its integrals over all lines in the 
plane (in two dimensions), or its integrals over all planes in space (in three dimensions). 
This has the consequence that in two-dimensional ACPAR, where the accessible data 
would in principle be the integrals over all lines in space (a four-dimensional manifold), 
typically only projections within a certain plane of rotation are considered, which reduces 
the three-dimensional problem to independent two-dimensional problems that can be 
solved by direct algorithms for the two-dimensional inverse Radon transform. While in 
the case of medical imaging such a sectioning approach is appropriate for minimizing 
the necessary radiation doses due to the elongated shape of the human body, this 
does not apply for projections of momentum densities in reciprocal space and thereby 
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corresponds to neglecting potentially independent information. In fact, due to the 
point group symmetry of the crystal taking only projections within a common rotation 
plane already implies information about out-of-plane projections in many cases, which a 
pure sectioning approach does not take into account. In a recent proposal (Kontrym- 
Sznajd & Samsel-Czekala 2007) this deficiency is addressed by re-parametrizing the 
reconstructed slices in terms of basis functions with the appropriate symmetry, but it is 
not known to which extent this can recover the information that was initially available, 
considering the correct symmetry, or whether the filtered reconstruction still reproduces 
the measured projections. In addition, the direct approaches typically need evenly and 
densely spaced projections in the plane of equal statistical precision, where the angular 
increment corresponds to the smallest resolvable features. In contrast, our proposed 
formulation allows an arbitrary number of projections with arbitrary orientation and 
noise level to be used, and takes into account all available information considering the 
prescribed symmetry. An early suggestion by Pecora (1987) to solve directly for the 
coefficients of three-dimensional basis functions of appropriate symmetry, dealing with 
the above-mentioned dimensional problem in a least-squares sense, should in principle 
be of equivalent quality in terms of these criteria, but apparently has seen only limited 
use. Note that it was shown there that it is not merely a welcome option to be able to 
use arbitrary orientations, but that using low-symmetry projections is actually indicated 
for optimal results in reconstruction. 

Direct transform methods are derived from analytical inversions of the mathematical 
Radon transform. As such, they cannot cope with experimental subtleties such as 
finite resolutions. Therefore, the resolution is typically deconvolved from the measured 
spectra before reconstruction. In most cases, maximum-entropy regularization is used 
to solve this ill-posed problem, either explicitly (Fretwell et al. 1995) or implicitly 
(Gerhardt et al. 1998), and sometimes even more arbitrary methods are used (Chiba 
et al. 2007). In contrast, aggregating all experimental complications (e.g., resolution, 
momentum sampling, shifting) into the forward operator as proposed here allows us 
to solve the inverse problem in a single step, subject to prior assumptions (equivalent 
to regularization) on the fundamental physical quantities, and thereby to rule out a 
compounding of the potentially conflicting regularization biases introduced at sequential 
steps. A very welcome additional consequence of our not touching the experimental 
spectra at all is that the noise statistics of the spectra remain unadulterated, specifically 
the Poissonian behaviour with independence between neighbouring pixels is conserved, 
which allows us to propagate the uncertainty of the data consistently into an estimated 
error of the final reported quantities. 

The effect of experimental noise on the reconstructed densities is a critical issue for 
the conventional methods. Especially with the direct transform methods that rely on the 
central slice theorem and interpolation (Kondo et al. 1993), the reconstruction is typically 
very ill-defined around the origin due to the conflicting information, and there is no 
obvious way of countering this problem. Such an effect is present also for basis expansion 
methods (Pecora 1987). In general, the regularity of the reconstructions is controlled by 
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the number of considered basis functions in the case of expansion methods, and by the use 
of an appropriate filter function in the various filtered transform methods, which in either 
case is rather opaque to the user. In contrast, in our formulation regularization by way 
of the prior distribution is explicit, and no noise artefacts appear in the reconstruction. 
Note that an explicit regularization functional is also the only practicable way to have 
non-linear biases such as a non-negativity constraint (as already observed by Pylak 
et al. 2011). 

In most published works, the identification of the Fermi surface is done subsequent 
to and independent from the reconstruction of the density. The simplest option is to 
transform the density from p- to k-space, i.e., subject it to the so-called LCW folding 
(Lock et al. 1973), and define the Fermi surface as an iso-density contour, e.g., according 
to a maximum gradient criterion (Biasini et al. 2002). The problem with such an 
approach in positron annihilation experiments is that due to the inequivalence between 
the measured two-photon momentum density and the electron momentum density, filled 
bands give rise to a non-constant background (Lock & West 1975), so that for finite 
resolution the Fermi surface does not strictly correspond to any iso-density contour. 
Edge-detection or enhancement methods (e.g., Dugdale et al. 1994, O’Brien et al. 1995) 
should be able to obviate this issue, but also in this case finite experimental resolution 
will tend to smooth regions of the Fermi surface with high curvature. Due to the explicit 
treatment of the Fermi surface during reconstruction, our proposed method does not- 
suffer from the above-mentioned effects. We are aware of only two comparable proposals: 
Biasini (2000) determined the parameters of a model Fermi surface by minimizing the 
deviation from the measured spectra after LCW folding, which is essentially the same idea 
as our proposal, only with constant band densities in k-space. This idea of reconstruction 
by a piecewise-constant function with explicit treatment of the step manifold has also 
been suggested in the context of diagnostical imaging (Ramlau & Ring 2007). Second, 
Laverock et al. (2010) propose to fit the LCW-folded spectra by calculated band densities, 
where the free parameters correspond to a state-dependent annihilation enhancement and 
energy shifts of the rigid bands. Clearly, in such an approach the freedom in the shape 
of the Fermi surface is very restricted, and the results will depend on the correctness of 
the electronic structure used as input. 

The only disadvantages of our proposed approach that we are aware of concern 
the increased numerical effort compared to direct methods, although this has ceased to 
be relevant with today’s computing power as we will show below. Also, no packaged 
software exists yet, but actually for momentum density reconstructions, different from 
medical imaging, typically custom implementations are used in any case, which especially 
in the case of orthogonal expansions can become more involved than our implementation 
as sketched below. 
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3. Implementation 

Due to the large dimensionality of the problem defined by Q, an efficient way to arrive 
at its solution is imperative for its practicability. We will give a detailed discussion of 
our implementation below, for its numerical aspects see the Supplementary Material. 

3.1. Parametrization of the solution space 

In our formulation, both the band densities p and the Fermi surface sheets cr are explicit 
degrees of freedom. As mentioned above, the physical picture suggests to describe 
the Fermi surface as the level sets of auxiliary smooth functions with the appropriate 
symmetry. For this purpose we use a Fourier description, where the required reciprocal- 
space translation symmetry is enforced by considering only those Fourier coefficients that 
correspond to real-space lattice vectors. Note that this is formally identical to a tight- 
binding description, and in fact it has been shown that it can reproduce the experimental 
Fermi surfaces in the noble metals with only a few free parameters (Roaf 1962). In 
contrast, for systems that conform rather to the free-electron picture, such as Al, different 
models will probably be more efficient in describing the Fermi surfaces (Ashcroft 1963). 

Thanks to the explicit treatment of the Fermi surface, the band densities p will 
be smooth and can therefore be described with comparatively low resolution. Here we 
use quadratic B-splines. The point symmetry is fulfilled by considering only coefficients 
within the irreducible wedge, with appropriate continuation over all of the reciprocal 
space. For reasons of efficiency we employ a non-constant sampling density, with higher 
resolution at low momentum. This is justified by the observation that the contribution 
at high momenta is mainly due to the core electrons, which will not be influenced much 
by the crystal structure and therefore be nearly isotropic. 

We want to emphasize that the choice of basis functions is not essential to the 
idea of the method. Considerable mathematical effort has been expended in deriving 
expansions that fulfill orthogonality relations of some sense under idealized projection 
operations (Louis 1984). In principle, such parametrizations could be used also in our 
proposal, where we expect that a sufficient description could be attained already at lower 
degrees of the expansion compared to previous implementations, as the parametrization 
needs to capture only the variation of the band densities excluding the Fermi surface 
steps. As an additional advantage, the system of linear equations to be solved would be 
better conditioned due to the orthogonality properties. However, such basis functions 
are typically non-vanishing nearly everywhere in reciprocal space, corresponding to full 
instead of sparse projection operators, which would constitute a serious drawback as we 
will discuss below. 

3.2. Matrix form of the operators 

The projection operator P specifies how a given three-dimensional momentum density 
leads to the set of two- or one-dimensional spectra in the specified orientations by line 
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or plane projections. For the construction of this matrix we compute for each voxel 
the set of pixels or bins in each spectrum that can potentially have an overlap with the 
projection of the voxel, and set the respective entry to the proportion received by the 
corresponding pixel. We find this proportion by an approach based on a look-up table 
with linear interpolation, pregenerated from a high-resolution projection of a single voxel 
cube. Given the finite experimental resolution (discussed below), this approach can be 
considered as equivalent to the exact solution. 

The optimal choice for the discretization of the three-dimensional momentum density 
is somewhat finer than the experimental resolution: A coarser representation will lead to 
artefacts, while a finer discretization will be numerically more expensive without any 
effect on the spectra after the application of resolution smearing. The same reasoning 
holds for the discretization of the spectra. As a consequence, each voxel will contribute 
to only a few pixels in each spectrum, allowing us to precompute and store the projection 
matrix in sparse format. 

By setting the resolution kernel to zero once its value has dropped below some 
threshold, also the resolution smearing operator has the form of a sparse matrix. The 
final step in obtaining the spectra from a given momentum density is the pointwise 
multiplication by the momentum sampling function, which can corresponds to a diagonal 
matrix. In principle, the projection operator P is then the product of these three 
precomputed sparse matrices, although due to the involved dimensionalities it is more 
efficient to utilize associativity of matrix multiplication in the further steps and never 
compute the product of the resolution smearing matrix and the projection matrix proper 
explicitly. 

The occupation operator X CT , computing the momentum density from the band 
densities and the Fermi surface a, is for fixed a represented as the horizontal concatenation 
of diagonal matrices, with the entries corresponding to the respective occupations. Here 
special care has to be taken with the voxels at the Fermi surface: For a continuous 
variation of the resulting spectra with the shape of the Fermi surface, the occupations 
have to be computed according to the proportions of the voxel within the Fermi surface 
(at the relevant temperatures the width of the Fermi surface compared to the voxel 
size is negligible). We accomplish this by the Gilat-Raubenheimer method (Gilat & 
Raubcnhcimer 1966), i.e., by linearizing the variation of the band energy within the 
voxel and computing the enclosed volume explicitly. 

Summing up, the linear operator A a that relates the band density parameter vector 
p to the spectra for a given shape of the Fermi surface cr is the matrix product PX CT of 
the matrices discussed above, modified by an additional degree of freedom corresponding 
to a constant background intensity. 

3.3. Algebraic solution of the Bayesian problem 

The fact that allows for an efficient computation of Q lies in the observation that 
puk(p,cr\y) can be approximated very accurately by a Gaussian distribution for fixed a. 
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Specifically, by equating the measured experimental spectra with their true value (i.e., 
the expected value before quantization) and expanding the logarithm of the Poissonian 
distribution around the maximum we have 

log(piik(p,<%)) = -§(A a p - y) T W(A a p -y) + const, (6) 

where the weighting matrix W is the inverse of the covariance matrix. This initial 
problem can be iteratively refined by expanding around the updated expected values. In 
our tests, such an iterative approach converged rapidly (within one iteration), and even 
the initial distribution was a faithful representation of the final result, as the counts per 
pixel were not too low and the residuals were small. 

To guarantee smooth reconstructed band densities, we formulate our prior 
distribution as a Gaussian distribution of the square norm of the band densities’ 
second derivatives plus an analogous contribution from the first derivatives to favour 
monotonicity. A non-constant weighting of these norms can be used to penalize the same 
density curvatures more if they occur at high momenta, where the absolute values of the 
densities are smaller. Our choice is functionally equivalent to Tikhonov regularization 
(Tikhonov 1963), i.e., a positive semidehnite quadratic form of the parameters, when 
searching for the maximum a posteriori estimate. 

As both pi ik and p pr i or are multivariate Gaussian distributions, so is their product 
Ppost- As a consequence, the maximum a posteriori estimate p*(cr) for given cr can be 
obtained by solving the system of linear equations 


(aJwa ct + £ a,d7 d Op = aJw y , 

i 


( 7 ) 


where Dj are the matrices that compute the (optionally weighted) derivatives of the 
band densities from the parameters and A; the corresponding regularization parameters. 
Also the marginal posterior distribution with respect to cr as defined by ([5]) follows easily 
as 


PpostHy) =P P ost(p*(o-),cr|y) 


/dettAjWA^ + ^AiDTDO. 


( 8 ) 


4. Application to model data 

Most previous proposals for algorithms to reconstruct three-dimensional momentum 
densities have been validated only by comparing the results on experimental data to 
those of other algorithms, if at all (e.g., Kontrym-Sznajd & Samsel-Czekala 2008, Pylak 
et al. 2011). Obviously this is not satisfactory, as specific aspects of the data could have 
been missed simultaneously by both the tested and the benchmark methods. In other 
cases the algorithms have been applied to synthetic data obtained by comparatively 
simplistic models, but also there the comparison has been done only in qualitative 
terms (e.g., Pecora 1987). As we claim here to be able to reconstruct Fermi surfaces 
quantitatively, we have to substantiate this claim by demonstrating both its fidelity (the 
magnitude of introduced systematic errors) and statistical performance (the propagation 
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Pon (r.l.u.) p 100 (r.l.u.) 

Figure 1. Simulated spectra for (100) (left) and (Oil) (right) orientation, including 
projection of the first Brillouin zone. The assumed resolution function is indicated. 


of experimental noise to the resulting dimensions) on realistic data for which the correct 
solution is known. 

For this purpose we chose the system of copper with its prototypical and well- 
known Fermi surface. We computed the electronic structure of Cu self-consistently in 
the generalized gradient approximation with the PBE exchange-correlation functional 
(Perdew et al. 1996) by the density functional code abinit (Gonze et al. 2009). With 
the converged density we computed both the electron wave functions and energies on a 
fine mesh and the T-point positron wave function (Barbiellini et al. 1995). Thanks to 
the plane-wave formulation used in the ABINIT code, the electron-positron momentum 
densities could be conveniently derived in a custom implementation corresponding to 
the independent particle model. From the three-dimensional density plus some constant 
background we computed the corresponding ACPAR spectra of 144 2 pixels for (001), 
(110) and (111) orientation, each with 25-10 6 counts distributed according to Poissonian 
statistics. We chose a discretization of 24 pixels per reciprocal lattice constant and 
assumed an anisotropic Gaussian resolution function with 2x1 pixel standard deviation. 
Two of the resulting spectra are given in Fig. [l] With the actual lattice constant of Cu, 
our chosen resolution would correspond to 1.32 x 0.66 mrad 2 FWHM. 

We reconstructed the density in a volume of 144 3 voxels at the discretization of 24 
voxels per reciprocal lattice constant. We want to emphasize that the coincidence of 
pixel and voxel size is by no means necessary; specifically for an actual experiment it 
will be beneficial to choose the voxel size as an integer fraction of the reciprocal lattice 
constant, while the discretization of the spectra will be given by the apparatus. As 
Cu has only a single band crossing the Fermi energy, we considered a fully occupied 
core level electron density and a single conduction band. We parametrized the space of 
band densities by cubic B-splines with an increasing sampling density towards small p, 
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p (r.l.u.) 


Figure 2. Selected cuts through the spectra: in (lOO)-spectrum along (Oil) (a), in 
(Oll)-spectrum along (Oil) (b) and (100) (c), in (lll)-spectrum along (Oil) (d) and 
(211) (e), in each case through the origin. 


corresponding to 385 degrees of freedom per band, and we described the single Fermi 
surface sheet by a five-parameter Fourier expansion (where the (110) coefficient is fixed 
to 1 as it corresponds to a trivial scaling of the energy range, and the (000) coefficient is 
chosen relative to the Fermi energy so as to constrain the occupied volume to half the 
Brillouin zone). We also added a term to the prior distribution to favour the decay of 
the coefficients with interaction range, because, as already noted by Roaf (1962), the 
Fermi surface changes only by very small amounts under certain modifications of the 
parameters. 

Maximizing the posterior probability given by Eq. ([8]) for the simulated spectra 
gives reconstructed spectra that, apart from the missing noise, are visually identical 
to the input spectra. This is substantiated by the one-dimensional cuts through the 
spectra shown in Fig. [ 2 ] and the corresponding reduced y 2 value of 1.047. The degree of 
achieved faithfulness to the data obviously depends on the choice of the regularization 
parameters X t in Eq. ([7]). Here we used the smallest values that still suppress visually 
noticeable artefacts (that would correspond to the reconstruction of experimental noise) 
in the reconstructed densities. 

A comparison of input and reconstructed density is given in Fig. [3j Here significant 
systematic differences can be discerned. Specifically, the original density is virtually 
constant within the Fermi surface in the first Brillouin zone and decays quite abruptly 
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Figure 3. Cuts through the three-dimensional electron-positron momentum density. 
Input densities due to model (left column) and reconstructed densities (right column), 
(100) (top row) and (Oil) (bottom row) planes through the origin, including outline of 
first Brillouin zone. 


through the necks into the second Brillouin zone. In contrast, the tendency of the 
reconstruction towards smooth variations in the band densities leads to a decrease 
already within the Fermi surface in the first Brillouin zone and to higher contributions 
from outer zones. Apart from that, the main features are clearly reconstructed in a 
qualitatively correct way. 

For assessing the algorithm’s performance in determining the Fermi surface in 
quantitative terms, we focus on three specific features: the extent of the Fermi surface 
along (100), along (110), and the radius of the (lll)-neck (its deviation from a circular 
shape is negligible). In the input data, these three dimensions are 1.063, 0.946 and 
0.203, respectively, measured in units of rf, the radius of the free-electron Fermi sphere. 
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In the reconstructions, the corresponding values are 1.063(3), 0.949(1) and 0.182(6), 
denoting the average and the standard deviations of the maximum-a-posteriori values for 
different realizations of the counting noise. For a single realization of the counting noise, 
the posterior probability distribution defined by Eq. ([8]) is to a good approximation a 
Gaussian distribution, with a covariance matrix that essentially corresponds to above- 
quoted standard deviations, which would allow us to estimate the errors of the dimensions 
obtained from an experiment. 

The small statistical uncertainty of the results indicates that our approach can 
also be used for qualitative statements. For instance, constraining the Fermi surface 
parameters in the reconstruction so that the necks along (111) become disconnected 
gives a maximum-a-posteriori probability that is smaller by 140 orders of magnitude, 
i.e., in a fictitious experiment this scenario could essentially be ruled out. 

5. Discussion 

The fact that our reconstruction reproduces the measured spectra perfectly within the 
errors shows that it considers all the information present in the data. Within these 
boundary conditions it yields the reconstruction most probable according to the prior 
assumptions. We think that our assumptions are arguably the soundest on physical 
grounds, and definitely the most transparent and easiest to adjust compared to those 
inherent to direct methods. 

Our results show that three ACPAR spectra with moderate statistics and 
unexceptional experimental resolution suffice for our interpretation method to give 
statistically very well-defined results on the Fermi surface dimensions. In terms of 
accuracy it has to be noted that specifically the radius of the (lll)-necks is underestimated 
by 10%. However, in the overall picture this corresponds only to an error of 2% of the 
mean Fermi surface radius, and all other regions of the Fermi surface are determined still 
much more accurate (see the juxtaposition of renderings of the input and reconstructed 
Fermi surface in Fig. [4]). For comparison, the Fermi surfaces of Ag and Au differ much 
more from the actual Fermi surface of Cu than the reconstruction does with a (lll)-neck 
of 0.137 in Ag and a (lOO)-radius of 1.135 in Au (Roaf 1962). Also, if a plausible model for 
the electronic structure in a given system is available, an analysis as presented here can be 
done on the model and the experimental results can be corrected in first approximation 
for the systematic effects displayed by the model reconstructions. Note further that 
also in quantum-oscillatory methods the Fermi surface has to be reconstructed from the 
measured data, subject to some plausibility assumptions. 

In the interpretation of an actual experiment it would probably be too optimistic 
to expect a ;\ /2 as low as reported here due to experimental imperfections such as 
additional contributions to the spectra, e.g., from surface positronium ejection, or slight 
misorientations. However, due to the adaptibility of our formulation such effects can be 
included at the cost of a few additional free parameters. This is in contrast to direct 
methods with, e.g., the strict assumption of crystal symmetry in the plane of sample 
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Figure 4. Rendering of the Fermi surface within first Brillouin zone, input data (left) 
and representative reconstruction from simulated spectra (right). 


rotation. 

After the alkali metals, the noble metals with Cu as the example chosen here have 
the simplest Fermi surface and are therefore probably the easiest systems to consider. 
For systems with multiple Fermi surface sheets separated only by a small distance the 
problem will become ill-defined due to resolution effects. Due to an analogous reasoning 
also large real space cells and consequently small reciprocal space features make a system 
hard to solve. However, these limitations obviously apply equally to any method of data 
interpretation. 

The last point we want to stress is that in principle any additional information 
can be considered in the prior probability assumptions. For example, in the case at 
hand the ill-defined (lll)-neck radius could be fixed to the value given by de Haas-van 
Alphen-measurements, where it corresponds to a prominent and well-defined frequency 
(Shoenberg 1962), while the shape of the Fermi surface belly is more directly encoded in 
the positron annihilation data. 

6. Conclusion 

Here we have presented a new point of view on obtaining the shape of the Fermi surface 
from Angular Correlation of Positron Annihilation Radiation, or equivalently Compton 
Scattering, data. We have pointed out the advantages of a unified formulation as an 
inverse problem in the Bayesian setting instead of the conventional sequential approaches. 
We have shown that modest requirements on the data statistics lead to statistically 
well-defined results on the Fermi surface dimensions, and we have discussed the small 
systematic biases introduced by the prior assumptions for the example of copper. These 
insights open the way for obtaining Fermi surface information with a quality that is 
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comparable to quantum-oscillatory methods under conditions such as high temperature 
and occupational disorder, where the classical methods cannot be applied. 
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